#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TStyle.h>
#include <TH1.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TF1.h>
#include <TTree.h>
#include <TChain.h>
#include <TFileCollection.h>
#include <TGraph2D.h>
#include <TLatex.h>
#include <TMarker.h>
#include <math.h>
#include <TMath.h>
#include <TFileInfo.h>

///Igor's stuff
#include "TROOT.h"
#include "MVAComputer2.h"
#include "PileUpProducer.h"
#include <TLorentzVector.h>
#include "MathRoot.h"
#include "InputVariablesHistograms.h"
//#define TESTIGORCODE


//#include "../interface/tiltt.h"

//#define SUSY
//#define MC
#define DATA


#undef PTMAX300

#ifdef MC
#define LUMIWEIGHTING
//#define HLTWEIGHTING
//#define bPARTONMATCH
#endif

#ifdef DATA
//#define MATCHING
#endif

#include "readhbb.h"

TCanvas* canvas;

const int nflav=3;
const int novl = 2;

class bTagEff {
public:
  std::string file;
  std::string onloff;
//   const int& nflav;
//   const int& novl;
  const std::string* flavlabel;
  const std::string* sovl;

  TFile* hfile;
  TH2D* hrel2D[nflav][novl];
  //TH2D*** hrel2D;


  bTagEff(){};
  bTagEff(const std::string& theFile,const std::string& theonloff,const std::string* theFlavlabel,
	  const std::string* theOvlLabel ){
    if (theonloff != "online" && theonloff != "offline") {
      std::cout << "bTagEff constructor: Bad online/offline mode " << theonloff << std::endl;
      return;
    }
    hfile = TFile::Open( theFile.c_str() );
    if (hfile == NULL) {
      std::cout << "bTagEff constructor: Cannot open file " << theFile << std::endl;
      return;
    }

    flavlabel = theFlavlabel;
    sovl = theOvlLabel;

    for (int iflav=0; iflav<nflav; ++iflav) {
      for (int iovl=0; iovl<novl; ++iovl) {
	if (theonloff == "online") {
	  hrel2D[iflav][iovl] = (TH2D*) hfile->Get( Form("pteta_%s_%s_releff",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
	  if (hrel2D[iflav][iovl] == NULL) std::cout << "Histogram not found: " 
			    << Form("pteta_%s_%s_releff",sovl[iovl].c_str(),flavlabel[iflav].c_str()) << std::endl;
	} else {
	  hrel2D[iflav][iovl] = (TH2D*) hfile->Get( Form("pteta_%s_%s_eff_div",sovl[iovl].c_str(),flavlabel[iflav].c_str()) );
	  if (hrel2D[iflav][iovl] == NULL) std::cout << "Histogram not found: " 
						     << Form("pteta_%s_%s_eff_div",sovl[iovl].c_str(),flavlabel[iflav].c_str()) << std::endl;
	}
      }
    }
  }
  double eff(int iflav,const std::string& theOvlLabel,double pt,double eta) {
    if (iflav<0 || iflav>=nflav) {
      std::cout << "eff: bad flavor " << iflav << std::endl;
      return -999999;
    }
    int iovl = -1;
    for (int jovl=0; jovl<novl; ++jovl) {
	if (sovl[jovl] == theOvlLabel) {
	  iovl = jovl;
	  break;
	}
    }
    if (iovl == -1) {
      std::cout << "Bad btag type requested " << theOvlLabel << std::endl;
      return -999999;
    }
#ifndef PTMAX300
    if (pt>300) pt=299.99;
#endif
    return hrel2D[iflav][iovl]->Interpolate(pt,fabs(eta));
  }
};





class TrigHistArray {
public:
  std::vector<TH1F*> histos;
  TH1F* histAllTrig;
  TH1F* histAllTrigWeighted;
  std::vector<std::string>* gtl;  // genericTriggerList
  unsigned int filterPattern;
  
  TrigHistArray() {
    std::cout << "Standard constructor called, this should not happen" << std::endl;
  }

  TrigHistArray(std::vector<std::string>* theGtl,std::vector<std::string>* theFilter,const char* genName,const char* genTitle,int nx,float xmin,float xmax) {
    gtl = theGtl;
    char ctn[1000];
    char ctt[1000];
    // first set the filter pattern
    filterPattern = 0;
    if (theFilter != NULL) {
      for (std::vector<std::string>::iterator fit=theFilter->begin(); fit != theFilter->end(); ++fit) {
	std::vector<std::string>::iterator tSlotBtag = std::find(theGtl->begin(), theGtl->end(), *fit);
	if (tSlotBtag != theGtl->end()) {
	  filterPattern = filterPattern | (1<<(tSlotBtag - theGtl->begin()));
	} else {
	  std::cout << "Filter trigger " << *fit << " not found in any slot" << std::endl;
	}
      }
    } else {
      filterPattern = ~0;  // set all filter bits to one
    }

    for (unsigned int ib=0; ib<gtl->size(); ++ib) {
      if (filterPattern & (1<<ib)) {
	sprintf(ctn,"%sTrig%d",genName,ib);
	sprintf(ctt,"%s Trigger %s",genTitle,(*gtl)[ib].c_str());
	TH1F* theHist = new TH1F(ctn,ctt,nx,xmin,xmax);
	if (theHist) {
	  histos.push_back(theHist);
	} else {
	  std::cout << "book; Problem creating histogram" << std::endl;
	}
      }
    }
    // inclusive histogram
    sprintf(ctn,"%sAllTrig",genName);
    sprintf(ctt,"%s all triggers",genTitle);
    histAllTrig = new TH1F(ctn,ctt,nx,xmin,xmax);
    sprintf(ctn,"%sAllTrigWeighted",genName);
    sprintf(ctt,"%s all triggers weighted",genTitle);
    histAllTrigWeighted = new TH1F(ctn,ctt,nx,xmin,xmax);
  }
			     
  void fill(int theTrgAccept,float x,float weight=1) {
    for (unsigned int ib=0; ib<gtl->size(); ++ib) {
      if (theTrgAccept & filterPattern & (1<<ib)) {
	histos[ib]->Fill(x,weight);
      }
    }
    histAllTrig->Fill(x);
    histAllTrigWeighted->Fill(x,weight);
  }
};


void getHbbMetadata(TFileCollection& fc,std::vector<std::string>& genericTriggerList,const float targetLumi,
		    float& lumiScaleFac,const bool doMC) {

  // extracts the following metadata from the Hbb analysis ntuples:
  //     list of generic triggers (= dictionary for bit number convention)
  //     total number of input events before any filtering
  //     generator cross sections and filter efficiency
  // computes lumi scaling factor corresponding to given target integrated luminosity

  lumiScaleFac = 1;
  int nInputEvents = 0;
  // first get the generic trigger list from the first file
  if (fc.GetNFiles() <=0) {
    std::cout << "getHbbMetadata: file list empty" << std::endl;
    return;
  }

  bool gtlDone = false;
  TIter iter( (TCollection*) fc.GetList() );
  while (TFileInfo* tfi = (TFileInfo*) &(*(iter.Next()))) {
    std::cout << "Open file " << tfi->GetCurrentUrl()->GetFile() << std::endl;
    std::cout << "Open Url " << tfi->GetCurrentUrl()->GetUrl() << std::endl;
    TFile* fg = TFile::Open( tfi->GetCurrentUrl()->GetUrl() );

    if (doMC) {
      TH1F* eventCount = (TH1F*) fg->Get("InputEvents/EventCount");
      if (eventCount != NULL) {
	std::cout << "Adding to event count " << eventCount->GetEntries() << std::endl;
	nInputEvents += eventCount->GetEntries();
      } else {
	std::cout << "eventCount histogram not found" << std::endl;
      }
    }

    if (! gtlDone) {
      TH1F* gtlHist = (TH1F*) fg->Get("hbbanalysis/gtlHist");
      if (gtlHist != NULL) {
	for (int ibin=1; ibin<= gtlHist->GetXaxis()->GetNbins(); ++ibin) {
	  genericTriggerList.push_back(gtlHist->GetXaxis()->GetBinLabel( ibin ));
	  std::cout << "getHbbMetadata: add generic trigger " << gtlHist->GetXaxis()->GetBinLabel( ibin )
		    << std::endl;
	}
	gtlDone = true;
      } else {
	std::cout << "gtlHist histogram not found" << std::endl;
      }
    }
    if (! doMC) break;  // for real data, only need first file to extract GTL
  }

  if (doMC) std::cout << "Total count of input events: " << nInputEvents << std::endl;
  std::cout << "Number of generic triggers read: " << genericTriggerList.size() << std::endl;

  if (doMC) {
    // now chain the GenInfo tree to access cross section & filter efficiency
    TChain gInfo("hbbanalysis/GenInfo");
    gInfo.AddFileInfoList((TCollection*) fc.GetList());
    TTree* gInfoTree = &gInfo;
    int nGInfo = gInfoTree->GetEntries();
    std::cout << "Number of events in gInfoTree: " << nGInfo << std::endl;

    // set GenInfo branches manually
    gInfoTree->SetBranchAddress("internalXsec",&_internalXsec);
    gInfoTree->SetBranchAddress("externalXsecLO",&_externalXsecLO);
    gInfoTree->SetBranchAddress("externalXsecNLO",&_externalXsecNLO);
    gInfoTree->SetBranchAddress("filterEfficiency",&_filterEfficiency);

    double SInternalXsec = 0;
    double SExternalXsecLO = 0;
    double SExternalXsecNLO = 0;
    double SFilterEfficiency = 0;
    for (int iGInfo=0; iGInfo<nGInfo; ++iGInfo) {
      gInfoTree->GetEntry(iGInfo);
      SInternalXsec += _internalXsec;
      SExternalXsecLO += _externalXsecLO;
      SExternalXsecNLO += _externalXsecNLO;
      SFilterEfficiency += _filterEfficiency;
    }
  
    double InternalXsec = SInternalXsec / nGInfo;
    double ExternalXsecLO = SExternalXsecLO / nGInfo;
    double ExternalXsecNLO = SExternalXsecNLO / nGInfo;
    double FilterEfficiency = SFilterEfficiency / nGInfo;

    std::cout << "Generator information:" << std::endl
	      << "      InternalXsec      = " << SInternalXsec / nGInfo << " pb " << std::endl
	      << "      ExternalXsecLO    = " << SExternalXsecLO / nGInfo << " pb " << std::endl
	      << "      ExternalXsecNLO   = " << SExternalXsecNLO / nGInfo << " pb " << std::endl
	      << "      FilterEfficiency  = " << SFilterEfficiency / nGInfo << std::endl;

    // in case not filter has been applied
    if (FilterEfficiency <0) {
      std::cout << "GenInfo FilterEfficieny <0 ==> set to 1" << std::endl;
      FilterEfficiency = 1;
    }

    // compute the lumi scale factor
    double XSec = InternalXsec;

///Igor's code
#ifdef SUSY
	XSec=111e0;
#endif

    double integLumi = nInputEvents / (XSec * FilterEfficiency);
    lumiScaleFac = targetLumi / integLumi;

    std::cout << "Lumi scale factor computation based on internal cross section:" << std::endl
	      << "      Xsec              = " << XSec << " pb " << std::endl
	      << "      IntLumiMC         = " << integLumi << " pb-1 " << std::endl
	      << "      TargetLumi        = " << targetLumi << " pb-1 " << std::endl
	      << "      Lumi Scale Factor = " << lumiScaleFac << std::endl;
  }
}

int jetFlavorCode(const int iJet) {
  int theFlavCode = -1;
  int theFlavor = partonFlavorJet[iJet];
  if (hflContentJet[iJet] != 0) theFlavor = hflContentJet[iJet];
  switch (abs(theFlavor)) {
  case 0:
  case 1:
  case 2:
  case 3:
    theFlavCode = 0;
    break;
  case 4:
    theFlavCode = 1;
    break;
  case 5:
    theFlavCode = 2;
    break;
  case 21:
    theFlavCode = 0;
    break;
  default:
    std::cout << "bad flavor " << theFlavor << std::endl;
  }
  return theFlavCode;
}


int diJetFlavorCode(const int iJet0, const int iJet1) {
  // determine dijet flavor from MC info
  // we set a flavor pair code that does not distinguish ordering
  int theFlav[2];
  theFlav[0] = jetFlavorCode( iJet0 );
  theFlav[1] = jetFlavorCode( iJet1 );
  // now set the dijet flavor
  int theFcDijet = -1;
  int flavPatternDijet = 10*theFlav[0]+theFlav[1];
  switch (flavPatternDijet) {
  case 0:
    theFcDijet = 5;
    break;
  case 11:
    theFcDijet = 3;
    break;
  case 22:
    theFcDijet = 0;
    break;
  case 1:
  case 10:
    theFcDijet = 4;
    break;
  case 2:
  case 20:
    theFcDijet = 2;
    break;
  case 12:
  case 21:
    theFcDijet = 1;
    break;
  default:
    std::cout << "diJetFlavorCode: Bad flavor code " << theFlav[0] << " " << theFlav[1] << std::endl;
  }
  return theFcDijet;
}


void readtribn(bool _useMVA=true) {

  // here we require three btagged jets


///Igor's code
Int_t _nevt_before_mva=0;
Int_t _nevt_after_mva=0;
Int_t _nevt_before_mva_trig=0;
Int_t _nevt_after_mva_trig=0;

Double_t _wgt_nevt=0;

  canvas = new TCanvas ("cg1","mycanvas",10,10,800,600);
  // open an ntuple file
  std::cout << " starting..." << std::endl;

  std::vector<std::string> genericTriggerList;
  std::vector<std::string> triggerFilterList;

#ifdef MC
  bool _doMC = true;
#else
  bool _doMC = false;
#endif

  // create the trigger filter list. Only histograms for these triggers will be created
  triggerFilterList.push_back("HLT_CentralJet46_BTagIP3D_CentralJet38_BTagIP3D");
  triggerFilterList.push_back("HLT_CentralJet46_CentralJet38_DiBTagIP3D_v1");

  // chain mode
  TChain f("hbbanalysis/HBBTo4B");
  TFileCollection fc("dum","","theMergeList.txt");

  // extract generic trigger list and number of input events from ntuple files
  float lumiScaleFac = 1;
  //getHbbMetadata(fc,genericTriggerList,1000.,lumiScaleFac,_doMC);
///Igor's code
getHbbMetadata(fc,genericTriggerList,501.,lumiScaleFac,_doMC);

  f.AddFileInfoList((TCollection*) fc.GetList());
  TTree* hbbtree = &f;

  // create the btag efficiency objects

  const std::string flavlabel[nflav] = {"udsg","c","b"};
  // define ovl labels
  const std::string sovl[novl] = { "TCHPT", "TCHP6" };
  // corresponding TCHP cuts
  double btcut[novl] = { 3.41, 6 };

  // define flavor pair classes
  const int nfcDijet = 6;
  string sfcDijet[nfcDijet] = { "bb", "bc", "bq", "cc", "cq", "qq" };
  const int nfc3rd = 3;
  string sfc3rd[nfc3rd] = { "q", "c", "b" };
  // flavor triplet classes
  const int nfcTrip = 6;
  string sfcTrip[nfcTrip] = { "bbb", "bbc", "bbq", "bcb", "bqb", "non-bb" };
  // the following are used for templates
  const int nfc = 3;  
  string sfc[nfc] = { "q", "c", "b" };
  const int ncateg = 3; // position of untagged jet
  const int ncorr = 2; // correction levels
  const int ntpat = 1; // trigger pattern (online btag), set to 8 later

//   bTagEff* bTagEffOffline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt-N14/pteta-q/plotmistag-b.root","offline",flavlabel,sovl);
  bTagEff* bTagEffOffline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt-N14/pteta-q-nb30/plotmistag-b.root","offline",flavlabel,sovl);



//   bTagEff* bTagReleffOnline = new bTagEff("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/results/job014/QCD_Pt_bEnriched-T14b/pteta-q/plotmistag-b.root","online",flavlabel,sovl);


  // open the bbPurity correction functions
  TF1* fbbfrac[novl][ncateg];
  TFile* bbPur = new TFile("/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/calib/bbPurity.root");
  for (int iovl=0; iovl<novl; ++iovl) {
    for (int icateg=0; icateg<ncateg; ++icateg) {
      fbbfrac[iovl][icateg] = (TF1*) bbPur->Get( Form("fbbfrac-ww-%s-Cat%d",sovl[iovl].c_str(),icateg) );
      if ( fbbfrac[iovl][icateg] == NULL ) {
	std::cout << "bbPur correctiobn function not found" << std::endl;
	return;
      }
    }
  }

  TFile* hout = new TFile("readtribn.root","recreate");

  //const bool oldAddDir = TH1::AddDirectoryStatus();
  TH1::AddDirectory(true);

  TH1::SetDefaultSumw2();

  Int_t nentries = (Int_t) hbbtree->GetEntries();
  cout << "Number of events in ntuple: " << nentries << endl;

#include "/afs/naf.desy.de/user/r/rmankel/scratch/HbbPat/CMSSW_4_2_3_patch2/src/DesyHiggsAnalyses/HBBAnalyses/test/readhbbset.C"


 // Histogram handler
  std::map<std::string, TH1F *> H1_;
  std::map<std::string, TH2F *> H2_;

  // counter histogram for flavor pair classes
  TH1F* hfc[novl];
  TH1F* hfcm[novl];
  TH1F* hfcww[novl];
  TH1F* hfcmww[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
    hfc[iovl] = new TH1F(Form("hfc%s",sovl[iovl].c_str()),Form("Flavor triple code 3* %s",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcm[iovl] = new TH1F(Form("hfcm%s",sovl[iovl].c_str()),Form("Flavor triple code 3* %s, mass 100-140",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcww[iovl] = new TH1F(Form("hfc%sww",sovl[iovl].c_str()),Form("Flavor triple code weighted as 3* %s",sovl[iovl].c_str()),19,-1.5,17.5);
    hfcmww[iovl] = new TH1F(Form("hfcm%sww",sovl[iovl].c_str()),Form("Flavor triple code weighted as 3* %s, mass 100-140",sovl[iovl].c_str()),19,-1.5,17.5);
    for (int ifcDijet=0; ifcDijet<nfcDijet; ++ifcDijet) {
      for (int ifc3rd=0; ifc3rd<nfc3rd; ++ifc3rd) {
	hfc[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcm[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcww[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
	hfcmww[iovl]->GetXaxis()->SetBinLabel(nfcDijet* ifc3rd + ifcDijet +2,Form("%s%s",sfcDijet[ifcDijet].c_str(),sfc3rd[ifc3rd].c_str()));
      }
    }
  }
  TH1F* hfctrip = new TH1F("hfctrip","Flavor condensed triple code",7,-1.5,6.5);
  for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
    hfctrip->GetXaxis()->SetBinLabel(ifcTrip+2,sfcTrip[ifcTrip].c_str());
  }

  const int minJet = 3;

  // book histogram arrays
  TrigHistArray* amdijetbt[novl];
  TrigHistArray* amdijetww[novl];
  TrigHistArray* aptjetbt[novl][minJet];
  TrigHistArray* aptjetww[novl][minJet];
  TrigHistArray* amdijetfcbt[novl][nfcTrip];
  TrigHistArray* amdijetfcww[novl][nfcTrip];
  TrigHistArray* aptjetfcbt[novl][nfcTrip][minJet];
  TrigHistArray* aptjetfcww[novl][nfcTrip][minJet];

  for (int iovl=0; iovl<novl; ++iovl) {
    amdijetbt[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
			     Form("mjjbt%s",sovl[iovl].c_str()),
			     Form("m(Jet1Jet2) 3*%s",sovl[iovl].c_str()),50,0,500 );
    for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
      amdijetfcbt[iovl][ifcTrip] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						     Form("mjjbt%sfc%s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
			     Form("m(Jet1Jet2) 3*%s fc %s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),50,0,500 );
    }
#ifdef MC
    amdijetww[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					Form("mjjww%s",sovl[iovl].c_str()),
					Form("m(Jet1Jet2) weighted as 3*%s",sovl[iovl].c_str()),50,0,500 );
    for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
      amdijetfcww[iovl][ifcTrip] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						     Form("mjjww%sfc%s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
			     Form("m(Jet1Jet2) weighted as 3*%s fc %s",sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),50,0,500 );
    }
#endif
    for (int iJ=0; iJ<minJet; ++iJ) {
      aptjetbt[iovl][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dbt%s",iJ,sovl[iovl].c_str()),
				   Form("pt leading jet %d 3*%s",iJ,sovl[iovl].c_str()),30,0,300);
      for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
	aptjetfcbt[iovl][ifcTrip][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dbt%sfc%s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
				   Form("pt leading jet %d 3*%s fc %s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),30,0,300);
      }
#ifdef MC
      aptjetww[iovl][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dww%s",iJ,sovl[iovl].c_str()),
				   Form("pt leading jet %d weighted as 3*%s",iJ,sovl[iovl].c_str()),30,0,300);
      for (int ifcTrip=0; ifcTrip<nfcTrip; ++ifcTrip) {
	aptjetfcww[iovl][ifcTrip][iJ] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				   Form("ptj%dww%sfc%s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),
				   Form("pt leading jet %d weighted as 3*%s fc %s",iJ,sovl[iovl].c_str(),sfcTrip[ifcTrip].c_str()),30,0,300);
      }
#endif
    }
  }

  TrigHistArray* anBtag = new TrigHistArray(&genericTriggerList,0,"nBtag","Number of online btags",10,-0.5,9.5 );
  TrigHistArray* aptJet[3];
  TrigHistArray* aetaJet[3];
  TrigHistArray* aphiJet[3];
  TrigHistArray* tcheJet[3];
  TrigHistArray* tchpJet[3];
  TrigHistArray* aIsJetWithBtag[3];

  for (int ii=0; ii<3; ++ii) {
    aptJet[ii] = new TrigHistArray(&genericTriggerList,0,Form("ptj%d",ii),Form("pt of %dth leading jet",ii+1),25,0,200);
    aetaJet[ii] = new TrigHistArray(&genericTriggerList,0,Form("aetaj%d",ii),Form("eta of %dth leading jet",ii+1), 100, -3.5, 3.5);
    aphiJet[ii] = new TrigHistArray(&genericTriggerList,0,Form("aphij%d",ii),Form("phi of %dth leading jet",ii+1), 100, -3.15, 3.15); 
    tcheJet[ii] = new TrigHistArray(&genericTriggerList,0,Form("tchej%d",ii),Form("TCHE of %dth leading jet, ptjet>50",ii+1), 60, -10, 20 );
    tchpJet[ii] = new TrigHistArray(&genericTriggerList,0,Form("tchpj%d",ii),Form("TCHP of %dth leading jet, ptjet>50",ii+1), 60, -10, 20 );
    aIsJetWithBtag[ii] = new TrigHistArray(&genericTriggerList,0,Form("isWithBtagj%d",ii),Form("IsJetWithBtag for %dth leading jet",ii+1),10,-3.5,6.5);
  }

  // templates
  TrigHistArray* amTemplate[nfc][novl][ncateg][ncorr][ntpat];
  for (int ifc=0; ifc<nfc; ++ifc) {
    for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	for (int icorr=0; icorr<ncorr; ++icorr) {
	  for (int itpat=0; itpat<ntpat; ++itpat) {
	    amTemplate[ifc][iovl][icateg][icorr][itpat]
	      = new TrigHistArray(&genericTriggerList,&triggerFilterList,
				  Form("mjjTemp_%s_%s_Cat%dCorr%dTpat%d",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat),
				  Form("mJet1Jet2 Template, %s %s Cat%d Corr%d Tpat%d",sfc[ifc].c_str(),sovl[iovl].c_str(),icateg,icorr,itpat),
				  50,0,500);
	  }
	}
      }
    }
  }
  
  // statistics for online btag pattern (double btag)
  TrigHistArray* atpat[novl][ncateg];
  TrigHistArray* atpatall[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	atpat[iovl][icateg] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
						Form("tpat_%s_Cat%d",sovl[iovl].c_str(),icateg),
						Form("Online btag trig pattern, 2*%s, Cat%d",sovl[iovl].c_str(),icateg),9,-0.5,8.5);
      }
      atpatall[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					 Form("tpat_%s_all",sovl[iovl].c_str()),
					 Form("Online btag trig pattern, 2*%s, all",sovl[iovl].c_str()),9,-0.5,8.5);
  }

  // statistics for online btag pattern (triple btag)
  TrigHistArray* atpattripall[novl];
  for (int iovl=0; iovl<novl; ++iovl) {
    atpattripall[iovl] = new TrigHistArray(&genericTriggerList,&triggerFilterList,
					   Form("tpattrip_%s_all",sovl[iovl].c_str()),
					   Form("Online btag trig pattern, 3*%s, all",sovl[iovl].c_str()),9,-0.5,8.5);
  }

#ifdef MC
  // double btag purity histograms
  TH1F* hmdibbt[nfcDijet][novl][ncateg];
  TH1F* hmdibww[nfcDijet][novl][ncateg];
  for (int ifcDijet=0; ifcDijet<nfcDijet; ++ifcDijet) {
    for (int iovl=0; iovl<novl; ++iovl) {
      for (int icateg=0; icateg<ncateg; ++icateg) {
	hmdibbt[ifcDijet][iovl][icateg] = new TH1F(Form("mdib_%s_%s_Cat%dbt",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
				   Form("Jet1Jet2 %s 2* %s Cat%d",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
						 50,0,500);
	hmdibww[ifcDijet][iovl][icateg] = new TH1F(Form("mdib_%s_%s_Cat%dww",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
				   Form("Jet1Jet2 %s weighted as 2* %s Cat%d",sfcDijet[ifcDijet].c_str(),sovl[iovl].c_str(),icateg),
						 50,0,500);
      }
    }
  }
#endif

  TrigHistArray* aMatchPattern = new TrigHistArray(&genericTriggerList,0,"matchPattern","matchPattern",20,-0.5,19.5);
  
  TrigHistArray* adphiJet = new TrigHistArray(&genericTriggerList,0,"adphij","delta-phi of two leading jets", 100, -3.15, 3.15); 
  
///Igor's code here
/// _MVA_ selection here: Ini of all need support
///It's done in the external driver program run.C
//gROOT->ProcessLine(".L MVACalculator.h+");
//gROOT->ProcessLine(".L MVAComputer.h+");
///Create MVAComputers


std::string mvaPath = "/afs/naf.desy.de/user/m/marfin/scratch/releases/CMSSW_4_2_3/src/AlpgenAnalysis/AlpgenAnalysis2/data/weights";
std::vector<std::string>  mvaMethods;
mvaMethods.push_back("Likelihood"); /// add any other methods if you are sure that they leave at mvaPath
MVAComputer2 * comp=0;
InputVariablesHistograms * MVAInput=0;
InputVariablesHistograms * MVAInputTrig=0;

if (_useMVA) comp = new MVAComputer2("MVAComputer2","MVAComputer2",mvaPath,mvaMethods);
if (_useMVA) MVAInput = new InputVariablesHistograms("Input","Input",mvaPath,mvaMethods);
if (_useMVA) MVAInputTrig = new InputVariablesHistograms("InputTrig","InputTrig",mvaPath,mvaMethods);

TH1 * _MVADiscr=0;
TH1 * _MVADiscrTrig=0;

if (_useMVA)  _MVADiscr = new TH1F("MVADiscr","MVADiscr",100,0.,1.);
if (_useMVA)  _MVADiscrTrig = new TH1F("MVADiscrTrig","MVADiscrTrig",100,0.,1.);

///_MVA_ ini stop



  // loop over the tree
  for (Int_t i=0; i<nentries /*&&  i<1e4  */; i++) {
    hbbtree->GetEntry(i);


if (i%1000==0) cout<<"Event # "<<i<<endl;
//cout<<"Event # "<<i<<endl;

    // this is mainly for counting events before any analysis cuts
    anBtag->fill(trgAccept,numberOfHltBtagObjects);
    
    double jetPtMin[3] = { 46, 38, 15};  // Roberval's cuts, but no extra cut on 3rd jet

#ifdef PTMAX300
    double jetPtMax[3] = {300, 300, 300}; // restrict to range where we have btag eff 
#else
    double jetPtMax[3] = {3500, 3500, 3500}; // restrict to range where we have btag eff 
#endif

    double maxEta = 2.6;

    std::vector<int> leadingJets;

///Igor's code
std::vector<math::XYZTLorentzVector> _vecJets ;
std::vector<double> * res;
double minPt=15e0; /// some cut on 

if (_useMVA) {

    // loop over the jets
     for (int iJet=0; iJet<numberOfJets; ++iJet) {
      if (! (fabs(etaJet[iJet])<maxEta) ) continue;	
      if (! (numberOfConstituentsInJet[iJet] > 1) ) continue;
      if ( ptJet[iJet]  < minPt ) continue;	


_vecJets.push_back(math::XYZTLorentzVector(pxJet[iJet], pyJet[iJet], pzJet[iJet], energyJet[iJet]));

	} /// loop over jets

if (comp) {
comp->ProcessEvent(_vecJets);
res = comp->GetOutPut();
} ///comp


if (MVAInput) MVAInput->FillHistograms(_vecJets);




if ( res && res->size()>0) 
if (_MVADiscr ) _MVADiscr->Fill(res->at(0));

///if (res->size()==0) cout<<"NO MVA OUTPUT!"<<endl;

} ///_useMVA

    int nJet = 0;
    // loop over the jets
    for (int iJet=0; iJet<numberOfJets; ++iJet) {
      if (nJet >= minJet) break;
      if (! (fabs(etaJet[iJet])<maxEta) ) continue;

      if (! (numberOfConstituentsInJet[iJet] > 1) ) continue;
      
      if ( (ptJet[iJet] > jetPtMin[nJet]) && (ptJet[iJet] < jetPtMax[nJet]) ) {
	leadingJets.push_back(iJet);
	++nJet;
      }
    }

    if (nJet < minJet) continue;

///End of preselection?




    // check the matching flags
    int mJet = 0;
    int matchPat = 0;
    for (int iJet=0; iJet<numberOfJets; ++iJet) {
      if (! (ptJet[iJet] > 15) ) continue;
      if (! (fabs(etaJet[iJet])<maxEta) ) continue;
      if (! (numberOfConstituentsInJet[iJet] > 1) ) continue;
      if (isJetWithHltBtag[iJet]) {
	matchPat = matchPat | (1<<mJet);
      }
      ++mJet;
    }
    aMatchPattern->fill(trgAccept,matchPat);

    // compute the weight for this event
    float weight = 1;
#ifdef LUMIWEIGHTING
    weight = lumiScaleFac;
#endif

    bool matchOK = true;
#ifdef MATCHING
    matchOK = false;
    if (nJet>=minJet ) {
      matchOK = isJetWithHltBtag[leadingJets[0]] && isJetWithHltBtag[leadingJets[1]];
    }
#endif

///Igor's code here

_nevt_before_mva++;
_wgt_nevt+=weight;

///Final preselection ?

bool _isOk=false;

float _deltaRj = -1;
    if ( (nJet>=minJet) && matchOK ) {

      float dphij =  phiJet[leadingJets[1]] - phiJet[leadingJets[0]];
      if (dphij>3.1415926) dphij -= 2*3.1415926;
      if (dphij<-3.1415926) dphij += 2*3.1415926;
      adphiJet->fill(trgAccept,dphij,weight);
      float detaj = etaJet[leadingJets[1]] - etaJet[leadingJets[0]];
      _deltaRj = sqrt( dphij*dphij + detaj*detaj );
    }

  bool _btag=false;
  if ( (tcHPBJetTag[leadingJets[1]]>btcut[0]) && (tcHPBJetTag[leadingJets[2]]>btcut[0]) ) {
   _btag = true;
  }
  if ((tcHPBJetTag[leadingJets[0]]>btcut[0]) && (tcHPBJetTag[leadingJets[2]]>btcut[0]) ) {
    _btag = true;
  }
  if ((tcHPBJetTag[leadingJets[0]]>btcut[0]) && (tcHPBJetTag[leadingJets[1]]>btcut[0]) ) {
    _btag = true;
  };


 if ( (nJet>=minJet) && (_deltaRj>1) && matchOK && _btag)  _isOk=true;

if  (_isOk)
if (trgAccept  & (1<<0) || trgAccept  & (1<<1) ) ///triggered data by Trigger0 or Trigger1
 _nevt_before_mva_trig++;

 

///_MVA_ selection here: Filtering
if (_useMVA) {
/// 1e-5 must to avoid artifact of the likelihood contruction, often  values ~10e-9
//double _cutMin=1e-8; /// No resrtriction on lower bound, mva output must be always positive 
double _cutMin=-1e0; /// No resrtriction on lower bound, mva output must be always positive

//double _cutMax=0.026; ///Tight cut
double _cutMax=0.18; ///Loose cut

/**
std::vector<math::XYZTLorentzVector> _vecJets ;
for (unsigned int k=0;k<leadingJets.size();k++) 
_vecJets.push_back(math::XYZTLorentzVector(pxJet[leadingJets[k]], pyJet[leadingJets[k]], pzJet[leadingJets[k]], energyJet[leadingJets[k]]));



std::vector<double> * res;
if (comp) {
comp->ProcessEvent(_vecJets);
res = comp->GetOutPut();
}
**/

#ifdef TESTIGORCODE
//std::cout<<"size of result "<<res->size()<<endl;
if (res && res->size()>0 ) std::cout<<"I'm event "<<i<<" with mva response "<<res->at(0)<<std::endl;
#endif


if ( res && res->size()>0) {

/**
if (_isOk)
if (_MVADiscr && res->at(0) >= _cutMin && res->at(0)<=1e0) _MVADiscr->Fill(res->at(0),weight);
**/

///Tiggered MVA output

if (_isOk)
if (trgAccept  & (1<<0) || trgAccept  & (1<<1) ) ///triggered data by Trigger0 or Trigger1
if (_MVADiscrTrig && res->at(0) >= _cutMin && res->at(0)<=1e0)  _MVADiscrTrig->Fill(res->at(0),weight);

//cout<<"_vecJets size "<<_vecJets.size()<<endl;
//cout<<"weight "<<weight<<endl;


/// Triggered input variables
if (_isOk)
if (trgAccept  & (1<<0) || trgAccept  & (1<<1) ) ///triggered data by Trigger0 or Trigger1
if (MVAInputTrig) MVAInputTrig->FillHistograms(_vecJets,weight);


///Main MVA selection
if (res->at(0) < _cutMin ||  res->at(0)>_cutMax)  continue;
}
#ifdef TESTIGORCODE
std::cout<<"I pass mva selection"<<std::endl;
#endif


}

if (_isOk)
if (trgAccept  & (1<<0) || trgAccept  & (1<<1) ) ///triggered data by Trigger0 or Trigger1
 _nevt_after_mva_trig++;

_nevt_after_mva++;

///_MVA_ selection stop






    // fill jet histograms
    float deltaRj = -1;
    if ( (nJet>=minJet) && matchOK ) {

      for (int iJ=0; iJ<minJet; ++iJ) {
	// jet kinematics
	aptJet[iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	aetaJet[iJ]->fill(trgAccept,etaJet[leadingJets[iJ]],weight);
	aphiJet[iJ]->fill(trgAccept,phiJet[leadingJets[iJ]],weight);
	aIsJetWithBtag[iJ]->fill(trgAccept,isJetWithHltBtag[leadingJets[iJ]],weight);
      }
      float dphij =  phiJet[leadingJets[1]] - phiJet[leadingJets[0]];
      if (dphij>3.1415926) dphij -= 2*3.1415926;
      if (dphij<-3.1415926) dphij += 2*3.1415926;
      adphiJet->fill(trgAccept,dphij,weight);
      float detaj = etaJet[leadingJets[1]] - etaJet[leadingJets[0]];
      deltaRj = sqrt( dphij*dphij + detaj*detaj );
    }

    // determine jet flavor from MC info
    int theFlav[minJet];
    for (int iJ=0; iJ<minJet; ++iJ) {
      theFlav[iJ] = jetFlavorCode(leadingJets[iJ]);
    }

    // determine the dijet flavor code for the first two leading jets
    int theFcDijet = diJetFlavorCode( leadingJets[0], leadingJets[1] );

    // now combine with flavor of third jet
    int theFc3rd = theFlav[2];
    int theFc = nfcDijet * theFc3rd + theFcDijet;

    // now make the triplet code
    int theFcTrip = -1;
    switch (theFcDijet) {
    case 0:  // (bb)
      switch (theFc3rd) {
      case 2: // bbb
	theFcTrip = 0;
	break;
      case 1: // bbc
	theFcTrip = 1;
	break;
      case 0: // bbq
	theFcTrip = 2;
	break;
      default:
	std::cout << "Bad 3rd flavor code " << theFcDijet << " " << theFc3rd << std::endl;
      }
      break;
    case 1:  // (bc)
      switch (theFc3rd) {
      case 2: // bcb
	theFcTrip = 3;
	break;
      case 1:
      case 0:
	theFcTrip = 5; // non-bb
	break;
      default:
	std::cout << "Bad 3rd flavor code " << theFcDijet << " " << theFc3rd << std::endl;
      }
      break;
    case 2:  // (bq)
      switch (theFc3rd) {
      case 2: // bqb
	theFcTrip = 4;
	break;
      case 1:
      case 0:
	theFcTrip = 5; // non-bb
	break;
      }
      break;
    case 3:
    case 4:
    case 5:
      theFcTrip = 5; // non-bb
      break;
    default:
      std::cout << "Bad flavor codes " << theFcDijet << " " << theFc3rd << std::endl;
    }

    hfctrip->Fill(theFcTrip,weight);

    // compute invariant mass of two leading jets
    if ( (nJet>=minJet) && (deltaRj>1) && matchOK ) {
      float energyTot = energyJet[leadingJets[0]] + energyJet[leadingJets[1]];
      float pxTot = pxJet[leadingJets[0]] + pxJet[leadingJets[1]];
      float pyTot = pyJet[leadingJets[0]] + pyJet[leadingJets[1]];
      float pzTot = pzJet[leadingJets[0]] + pzJet[leadingJets[1]];
      
      float dijet_mass_sq = energyTot * energyTot - pxTot * pxTot - pyTot * pyTot - pzTot * pzTot;

      if (dijet_mass_sq >= 0) {
	float dijet_mass = sqrt( dijet_mass_sq );

	// fill the triple-btag histograms
	for (int iovl=0; iovl<novl; ++iovl) {
	  if ( (tcHPBJetTag[leadingJets[0]]>btcut[iovl]) &&  (tcHPBJetTag[leadingJets[1]]>btcut[iovl])
	       && (tcHPBJetTag[leadingJets[2]]>btcut[iovl]) ) {
	    amdijetbt[iovl]->fill(trgAccept,dijet_mass,weight);
	    // triple-jet specific flavor code histograms
	    if (theFcTrip>=0) amdijetfcbt[iovl][theFcTrip]->fill(trgAccept,dijet_mass,weight);
	    // fc counting
	    hfc[iovl]->Fill( theFc,weight );
	    if ( (dijet_mass > 100) && (dijet_mass < 140 ) ) hfcm[iovl]->Fill( theFc,weight );
	    for (int iJ=0; iJ<minJet; ++iJ) {
	      aptjetbt[iovl][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	      if (theFcTrip>=0) aptjetfcbt[iovl][theFcTrip][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],weight);
	    }
	    // tpat counting
	    int theTpat = 4 * int(isJetWithHltBtag[leadingJets[2]]) + 2 * int(isJetWithHltBtag[leadingJets[1]])
	      + int(isJetWithHltBtag[leadingJets[0]]);
	    atpattripall[iovl]->fill(trgAccept,float(theTpat),weight);
	  }

#ifdef MC
	  // now do the same without btag cut and with weighting
	  double wbt[minJet];
	  for (int iJ=0; iJ<minJet; ++iJ) {
	    wbt[iJ] = bTagEffOffline->eff(theFlav[iJ],sovl[iovl].c_str(),ptJet[leadingJets[iJ]],etaJet[leadingJets[iJ]]);
	  }
	  double wtot = weight * wbt[0] * wbt[1] * wbt[2];
	  amdijetww[iovl]->fill(trgAccept,dijet_mass,wtot);
	  // triple-jet specific flavor code histograms
	  if (theFcTrip>=0) amdijetfcww[iovl][theFcTrip]->fill(trgAccept,dijet_mass,wtot);
	  // fc counting
	  hfcww[iovl]->Fill( theFc,wtot );
	  if ( (dijet_mass > 100) && (dijet_mass < 140 ) ) hfcmww[iovl]->Fill( theFc,wtot );
	  for (int iJ=0; iJ<minJet; ++iJ) {
	      aptjetww[iovl][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],wtot);
	      if (theFcTrip>=0) aptjetfcww[iovl][theFcTrip][iJ]->fill(trgAccept,ptJet[leadingJets[iJ]],wtot);
	  }
#endif
	}

	// background templates

	for (int iovl=0; iovl<novl; ++iovl) {
	  // set category membership for this btag selection
	  bool lcateg[3] = {false, false, false};
	  if ( (tcHPBJetTag[leadingJets[1]]>btcut[iovl]) && (tcHPBJetTag[leadingJets[2]]>btcut[iovl]) ) {
	    lcateg[0] = true;
	  }
	  if ((tcHPBJetTag[leadingJets[0]]>btcut[iovl]) && (tcHPBJetTag[leadingJets[2]]>btcut[iovl]) ) {
	    lcateg[1] = true;
	  }
	  if ((tcHPBJetTag[leadingJets[0]]>btcut[iovl]) && (tcHPBJetTag[leadingJets[1]]>btcut[iovl]) ) {
	    lcateg[2] = true;
	  };
	  // determine the online btag pattern
	  int theTpat = 4 * int(isJetWithHltBtag[leadingJets[2]]) + 2 * int(isJetWithHltBtag[leadingJets[1]])
	    + int(isJetWithHltBtag[leadingJets[0]]);
	  for (int icateg=0; icateg<ncateg; ++icateg) {
	    if (lcateg[icateg]) atpat[iovl][icateg]->fill(trgAccept,float(theTpat),weight);
	  }
	  atpatall[iovl]->fill(trgAccept,float(theTpat),weight);
// 	  if (lcateg[0] || lcateg[1] || lcateg[2]) {
// 	    std::cout << isJetWithHltBtag[leadingJets[0]] << " "
// 		      << isJetWithHltBtag[leadingJets[1]] << " "
// 		      << isJetWithHltBtag[leadingJets[2]] << " "
// 		      << "tPat = " << theTpat 
// 		      << "  categ " << lcateg[0] << lcateg[1] << lcateg[2] << std::endl;
// 	  }
	  // fill the templates
	  for (int ifc=0; ifc<nfc; ++ifc) {
	    for (int icateg=0; icateg<ncateg; ++icateg) {
	      if (lcateg[icateg]) {
		float mistagWeight = bTagEffOffline->eff(ifc,sovl[iovl].c_str(),ptJet[leadingJets[icateg]],etaJet[leadingJets[icateg]]);
		for (int icorr=0; icorr<ncorr; ++icorr) {
		  double wbbpur = 1;
		  if (icorr == 1) wbbpur =  fbbfrac[iovl][icateg]->Eval( dijet_mass );
		  for (int itpat=0; itpat<ntpat; ++itpat) {
		    // no additional corrections or trigger patterns yet
		    amTemplate[ifc][iovl][icateg][icorr][itpat]->fill(trgAccept,dijet_mass,weight * mistagWeight * wbbpur);
		  }
		}
	      }
	    }
	  }
	}

#ifdef MC
	// here we fill the double btag purity histograms
	// auxiliary array to perform permutations among first three jets
	int shift[6] = { 0, 1, 2, 0, 1, 2 };
	
	for (int iovl=0; iovl<novl; ++iovl) {
	  for (int icateg=0; icateg<ncateg; ++icateg) {
	    int jetA = leadingJets[shift[icateg+1]];
	    int jetB = leadingJets[shift[icateg+2]];
	    int dibFcDijet = diJetFlavorCode( jetA, jetB );
	    if ( (tcHPBJetTag[jetA]>btcut[iovl]) &&  (tcHPBJetTag[jetB]>btcut[iovl]) ) {
	      hmdibbt[dibFcDijet][iovl][icateg]->Fill( dijet_mass, weight );
	    }
	    // same with btag weighting
	    hmdibww[dibFcDijet][iovl][icateg]->Fill( dijet_mass, weight
						     * bTagEffOffline->eff(theFlav[shift[icateg+1]],sovl[iovl].c_str(),ptJet[jetA],etaJet[jetA])
						     * bTagEffOffline->eff(theFlav[shift[icateg+2]],sovl[iovl].c_str(),ptJet[jetB],etaJet[jetB]) );
	  }
	}

#endif
      } else {
	std::cout << "Unphysical dijet_mass_sq = " << dijet_mass_sq << std::endl;
	//dijet_mass = -1;
      } // if (dijet_mass...)
    } // if (njet...)
  }


///Igor's code
cout<<"MVA selection :: Total # of events  "<<nentries<<endl;
cout<<"MVA selection :: Total # of weighted events  "<<_wgt_nevt<<endl;
cout<<"MVA selection :: # Events come to MVA computer "<<_nevt_before_mva<<endl;
cout<<"MVA selection :: # Events come from MVA computer "<<_nevt_after_mva<<endl;
if (_nevt_after_mva>0) cout<<"MVA selection ::  Eff of  MVA computer "<<(Double_t)_nevt_after_mva/_nevt_before_mva<<endl;

cout<<"MVA selection :: # Triggered Events come to MVA computer "<<_nevt_before_mva_trig<<endl;
cout<<"MVA selection :: # Triggered Events come from MVA computer "<<_nevt_after_mva_trig<<endl;
if (_nevt_after_mva>0) cout<<"MVA selection ::  Eff of  MVA computer (Triggered Events) "<<(Double_t)_nevt_after_mva_trig/_nevt_before_mva_trig<<endl;



  // termination
  hout->Write();
  hout->Close();

if (MVAInputTrig) {
	delete MVAInputTrig;
}
if (MVAInput) {
delete MVAInput;
}


}
